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Abstract 

Today, the 'hydrogen atom model' is known to play its role not only in teaching 
the basic elements of quantum mechanics but also for building up effective theories in 
atomic and molecular physics, quantum optics, plasma physics, or even in the design of 
semiconductor devices. Therefore, the analytical as well as numerical solutions of the 
hydrogen-like ions are frequently required both, for analyzing experimental data and for 
carrying out quite advanced theoretical studies. In order to support a fast and consistent 
access to these (Coulomb-field) solutions, here we present the Dirac program which has 
been developed originally for studying the properties and dynamical behaviour of the 
(hydrogen-like) ions. In the present version, a set of Maple procedures is provided for 
the Coulomb wave and Green's functions by applying the (wave) equations from both, the 
nonrelativistic and relativistic theory. Apart from the interactive access to these functions, 
moreover, a number of radial integrals are also implemented in the Dirac program which 
may help the user to construct transition amplitudes and cross sections as they occur 
frequently in the theory of ion-atom and ion-photon collisions. 
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PROGRAM SUMMARY 



Title of program: DlRAC 
Catalogue number: To be assigned. 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland. 
Licensing provisions: None. 

Computer for which the program is designed and has been tested: All computers with a license 
of the computer algebra package Maple [1]. 

Program language used: Maple 8 and 9. 

Memory required to execute with typical data: 

No. of bytes in distributed program, including test data, etc.: 

Distribution format: tar gzip file 

CPC Program Library Subprograms required: None. 

Keywords: Analytical solution, Coulomb-Green's function, Coulomb problem, Dirac equation, 
energy level, expectation value, hydrogen-like ion, hydrogenic wavefunction, matrix element, 
radial integral, special functions. 

Nature of the physical problem: 

Analytical solutions of the hydrogen atom are widely used in very different fields of physics 
[2,3]. Despite of the rather simple structure of the hydrogen-like ions, however, the underlying 
'mathematics' is not always that easy to deal with. Apart from the well-known level structure 
of these ions as obtained from either the Schrddinger or Dirac equation, namely, a great 
deal of other properties are often needed. These properties are related to the interaction of 
bound electron(s) with external particles and fields and, hence, require to evaluate transition 
amplitudes, including wavefunctions and (transition) operators of quite different complexity. 
Although various special functions, such as the Laguerre polynomials, spherical harmonics, 
Whittaker functions, or the hypergeometric functions of various kinds can be used in most 
cases in order to express these amplitudes in a concise form, their derivation is time consuming 
and prone for making errors. In addition to their complexity, moreover, there exist a large 
number of mathematical relations among these functions which are difficult to remember in 
detail and which have often hampered quantitative studies in the past. 

Method of solution: 

A set of Maple procedures is developed which provides both the nonrelativistic and rela- 
tivistic (analytical) solutions of the 'hydrogen atom model' and which facilitates the symbolic 
evaluation of various transition amplitudes. 
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Restrictions onto the complexity of the problem: 

Over the past decades, a large number of representations have been worked out for the hydro- 
genic wave and Green's functions, using different variables and coordinates [2]. From these, 
the position-space representation in spherical coordinates is certainly of most practical inter- 
est and has been used as the basis of the present implementation. No attempt has been made 
by us so far to provide the wave and Green's functions also in momentum space, for which 
the relativistic momentum functions would have to be constructed numerically. 
Although the Dm AC program supports both symbolic and numerical computations, the latter 
one are based on Maple's standard software floating-point algorithms and on the (attempted) 
precision as defined by the global Digits variable. Although the default number, Digits = 
10, appears sufficient for many computations, it often leads to a rather dramatic loss in the 
accuracy of the relativistic wave functions and integrals, mainly owing to Maple's imprecise 
internal evaluation of the corresponding special functions. Therefore, in order to avoid such 
computational difficulties, the Digits variable is set to 20 whenever the Dirac program is 
(re-)loaded. 

Unusual features of the program: 

The Dirac program has been designed for interactive work which, apart from the standard 
solutions and integrals of the hydrogen atom, also support the use of (approximate) semirela- 
tivistic wave functions for both, the bound- and continuum states of the electron. To provide 
a fast and accurate access to a number of radial integrals which arise frequently in applica- 
tions, the analytical expressions for these integrals have been implemented for the one-particle 
operators r k , e~ ar , d m /dr m , jh{kr) as well as for the (so-called) two-particle Slater integrals 
which are needed to describe the Coulomb repulsion among the electrons. Further procedures 
of the Dirac program concern, for instance, the conversion of the physical results between 
different unit systems or for different sets of quantum numbers. A brief description of all pro- 
cedures as available in the present version of the Dirac program is given in the user manual 
Dirac-commands.pdf which is distributed together with the code. 

Typical running time: Although the program replies promptly on most requests, the running 
time also depends on the particular task. 

References: 

[1] Maple is a registered trademark of Waterloo Maple Inc. 

[2] H. A. Bethe and E. E. Salpeter, Quantum Mechanics of One- and Two-Electron Atoms, 
Springer, Berlin, 1957. 

[3] J. Eichler and W. Meyerhof, Relativistic Atomic Collisions, Academic Press, New York, 
1995. 
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LONG WRITE-UP 



1 Introduction 

Since the beginning of quantum mechanics, the hydrogen atom has been one of the best 
studied models in physics ^ n a iar ge number of textbooks and monographs, therefore, 

its (well-known) solutions from either the nonrelativistic Schrodinger equation or Dirac's 
(relativistic) wave equation have served as standard examples in introducing atomic theory. 
But apart from the educational insight as obtained from the hydrogen atom model, the theory 
of the one-electron atoms and ions has found its way also into quite different fields of modern 
physics including molecular and plasma physics, quantum optics, or the theory of quantum 
information. When combined with the atomic shell model, i.e. the successive filling of shells 
as established for instance by the periodic table of elements, the (analytical) solutions of the 
'hydrogen atom' help understand many atomic processes in Nature, at least qualitatively. 
During the last decades, moreover, Dirac's theory of the hydrogen-like ions has been applied 
successfully for describing energetic collisions of high-Z ions with atomic targets |3| or for 
testing quantum-electro dynamical (QED) effects [3], if Furry's picture and field-theoretical 
concepts are taken into account. 

As seen from this wide-spread use of the hydrogen atom model, therefore, a fast and interactive 
access to the solutions and properties of the one-electron ions is highly desirable. In the past, 
several program have been developed and implemented into the CPC program library for 
calculating the hydrogenic wave and Green's functions, both within the nonrelativistic El 
and relativistic jTllH] framework. Beside of providing these functions explicitly, they have been 
implemented also — at least to a certain extent — in a large number of further codes including 
those for studying the two-photon ionization of hydrogen-like ions jH] , their interaction with 
lasers |S] and electrons jlUj , the excitation and ionization of such ions by Coulomb fields jllj , 
their bound-bound transitions in the presence of external magnetic fields ^21) an d at many 
places elsewhere. But although these programs have certainly found useful for understanding 
the properties and behaviour of one-electron atoms and ions, they have two drawbacks in 
common which restrict their further application: Since these programs are designed for a very 
particular task, they often (i) cannot be extended so easily to other problems. Moreover, by 
using standard languages such as FORTRAN or C, these programs do (ii) not support any 
algebraic (symbolic) access to the wave functions and matrix element as required by more 
advanced studies. The lack of symbolic techniques, in particular, and together with the low 
flexibility in dealing with the input and output of the programs from above have hampered 
their use for a much wider range of applications. 

Today, a new and promising alternative for having an algebraic access to complex (symbolic) 
expressions is offered by computer algebra including the general purpose systems such as 
Mathematica or Maple. In the last decade, they have been increasingly utilized to provide 
a fast and reliable gateway to a large set of complex expressions and computations. With the 
development of the Dirac program, we here provide a set of Maple procedures for studying 
the properties and the dynamical behaviour of the hydrogen-like ions. In the past two years, 
this program has been used to support not only the interpretation of experiments but also 
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for a detailed analysis of the polarization and correlation properties of the x-ray radiation as 
emitted by relativistic collisions |131ll4j . Owing to its interactive design, moreover, the Dirac 
program may help the user not only in his or her daily research work but also in teaching the 
physics of (hydrogen-like) atoms and ions. 

Obviously, however, a program about the hydrogen atom model can neither be implemented 
comprehensively nor explained within a single step. Therefore, the Dirac package will be 
presented in several steps, each of them dedicated to a well defined part of the theory. To 
provide access to the very basic elements of the theory, here we first describe a set of proce- 
dures for the symbolic (as well as numerical) computation of the wave and Green's functions. 
These functions are provided for both, the bound- and continuum-state solutions of the non- 
relativistic as well as relativistic theory. Having once a simple access to these functions, of 
course, they can be rather easily applied to construct the transitions amplitudes and matrix 
elements for a large number of other properties. For several — frequently required — matrix 
elements, an algebraic evaluation of the radial integrals has been implemented also explicitly. 

Although the main purpose of this work is to describe the design and the use of the Dirac 
package, in the next section we first compile the basic equations and expressions from the 
theory of the hydrogen atom with emphasis on those which have been implemented in the code. 
In Section |31 we later describe the organization of the program and how it is distributed. This 
sections also lists all the user-accessible commands together with a short description, while a 
detailed explanation of the parameters, optional arguments and some additional information 
for each procedure is given in an additional user manual and is appended to the code. This is 
followed in Section 0] by a number of selected examples to illustrate the interactive use of the 
Dirac program within its Maple environment. A short outlook onto future developments 
and extensions of the program is finally given in Section \5\ 

2 Theoretical background 

Since the quantum-mechanical theory of the one-electron atoms has been worked out long 
time ago, both within the nonrelativistic and relativistic framework, we may refer for most 
details to the literature. Helpful introductions into this topic can be found, for instance, in 
the texts by Bethe and Salpeter PQ, Messiah |15| . Bransden and Joachain Eichler and 
Meyerhof 0|, and by many others. In the present work, therefore, we will restrict ourselves 
to rather a short compilation of the basic expressions, just enough in order to introduce the 
notations and formulas as implemented within the Dirac program. In the following two 
subsections 12.11 and 12.21 in particular, we briefly recall the (analytically well known) radial- 
angular representation of the hydrogen wave and Green's functions which are later applied in 
subsection 12.31 for evaluating a number of radial integrals. Emphasis has first been paid to 
those integrals which are of frequent use in the computation of the structural and dynamical 
properties of the hydrogen-like ions. In this paper (and also in most parts of the underlying 
implementation of the Dirac code), atomic Hartree units h = m e = e 2 /Aire® = 1, where eo 
is the permittivity of free space, are used throughout unless specified otherwise. 
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2.1 Wavefunctions and energy levels 
2.1.1 Nonrelativistic framework 

In a pure Coulomb potential V(r) = —Zjr of a given (nucleus with) charge Z, the motion of 
an electron is described by the time-independent Schrodinger equation 

A^(r) + 2 (e + ^) i>(r) = (1) 

which, in polar coordinates r, and <p, can be separated into three independent equations 
owing to the ansatz 

i/>(r,0,(p) = R(r)Y lm (d,<p) = ^-Y lm ($,ip) . (2) 

r 

As seen from this ansatz Q, the electronic wavefunction il>(r) = ip(r,6,(j)) is therefore given 
as a product of a radial function R(r) or P(r)/r, respectively, and a spherical harmonic 
Yi m ('d,ip) which describes the angular structure and which depends on the orbital angular 
momentum I and its projection onto the quantization axis as denoted by the magnetic quantum 
number m. However, not much need to be said here about the properties and treatment of 
the spherical harmonics as one of the standard tasks from the theory of angular momentum 
|16j for which meanwhile also a number of symbolic tools are available [171 118j . The (radial) 
function P(r) from ansatz © obeys the radial Schrodinger equation 



1 d f 2 d\ 1(1 + 1) | 2Z | ^ 



P(r) 

{) -0 (3) 



r 

and are known to belong either to the discrete part of the spectrum with the (negative) 
energies 



Z 



2 



En = -2^2 < ' n = 1 ' 2 '-' ( 4 ) 

the so-called bound states, or to its continuous part for all energies E > 0, for which the 
functions P(r) are said to represent the free-electron (or continuum) states. For a Coulomb 
potential in Eq. ©, an analytic solution of the radial functions P(r) in terms of the confluent 
hypergeometric function F(a,/3;x) are known for both, the bound-states 

,. , 1 (/' + ')! /-Z\ 3/2+ ' (_Zr) 

P m(.r) = r l+1 rni ; — K - ( — ) el J 

2Zr\ 

l-l), 21 + 2, , (5) 

n J 

with n and I being the principal and the orbital angular momentum quantum numbers, as 
well as for the continuum (E > 0) 

P El (r) = II ^ s2 + n ' 2 7§$T)T e ~ lkr F( - in ' + l + + ^ 2ifcr) (6) 

with k = \[2E and n' = Z/k. In the standard theory, moreover, the radial functions 
and © are normalized due to the conditions 

/>oo 

R 2 nl (r)r 2 dr = / P^r) dr = 1 (7) 
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for having a single particle per bound state and due to 

/•oo 

/ PUr)Pm{r)dr = 5{E-E') (8) 
J o 

for representing one particle per unit energy as far as electrons in the continuum are concerned. 

In the DlRAC package, both the full wavefunctions (j2J and their radial components © and © 
are provided for the bound- and free-electron states by using the procedures Dirac_orbital() 
and Dirac_radial_orbital(), respectively. 



2.1.2 Relativistic framework 

Similar to Schrodinger's equation the Dirac (eigenvalue) equation with a Coulomb poten- 
tial 

ica-V + - - (B - 1) c 2 + E J r/>(r) = (9) 
r / 

can be used to describe the relativistic motion of a spin-1/2 electron in the field of a nucleus 
with charge Z which, apart from the kinetic and potential energy of the electron, now also 
contains its rest energy as well as contributions due to the electron spin 012]. Analogue to 
the nonrelativistic case, the solution ip(r) of Eq. © can be separated again into a radial and 
angular part 



1 



P(r) n Km (#, <p) 



fl>Km{r,#,<p) = - , (10) 



r 



iQ(r) 0_ Km (#, <p) 



where f2 Km (i?, <p) denotes a standard Dirac spin-orbital function and where k = ± ( j + 1/2) 
for I = j ± 1/2 is called the relativistic angular momentum quantum number. Owing to the 
definition of k = ±1, ±2, this quantum number carries information about both, the total 
angular momentum j and the parity (-1)' of the wavefunction. As usual, the Dirac spin- 
orbitals can be written in terms of the spin-1/2 Pauli spinors x± an d t ne spherical harmonics 
by 

^m(^) = ^ {I, mi, 1/2, m - mi \ jm) Yi mi (ft, (p)xm- mi ■ (H) 
mi 

As seen from ansatz (|l(jp. moreover, the radial part of the wavefunctions is now given by the 
two functions P(r) and Q(r) which are called the large and small components, respectively. 
Similar to the nonrelativistic case, an analytic representation of these two components is 
known and is different for an electron in a bound or continuum state where, apart from 
the continuum with positive energies E > an additional continuous branch (of solutions) 
with negative energies E < —2c 2 exists. As known from Dirac's theory, these solutions with 
negative energies require special care and can be re-interpreted, for instance, in terms of 
positron states. For the bound states of hydrogen-like ions, the (discrete eigen-) energies of 
Eq. @ are given by 

F (W nK - 1) _ 1 J_ 

nK ~ ~ - , 2 -|V2 a 2 ' 1 ' 

aZ 



a 2 



1 + 



7 



while the corresponding radial components P nK and Q nK can be written again in terms of the 
confluent hyper geometric function F(a, (3; x). They read as 



Pn.(r) = N nK y/1 + W nK r (2qr) s - 1 e~ qr [-ri F{-n' + 1, 2s + 1; 2qr) 

aZ 



F(-ri,2s + l;2qr) 



(13) 



Q nK {r) = -N nK y/l - W nK r (2qr) s - 1 e^ r [ri F(-n' + l,2s + 1; 2qr) - 

olZ\ „, , 
r- F(-ri,2s + l;2qr) 

qKJ 



(14) 



where ri = n — \k\ = 0, 1, 2, ... denotes the number of nodes of the radial components, 
A c = h/m e c the Compton length of the electron, and 



s = \J k 2 — (aZ) 2 , 
Z 

9 ~ y/(aZ) 2 + (ri TW 
Moreover, the normalization factor 

V2q 5 / 2 X C 



(15) 
(16) 



r(2a + l) 



T(2s + ri + 1) 
n'\ (aZ) (aZ — nq\ c 



1/2 



(17) 



is chosen in such a way that the two radial components satisfy the condition 

/•oo 

/ (f«(r) + QL(r)) dr = 1, 
J o 



(18) 

i.e. that the wavefunction (|10|) represents one electron per bound state. 

For the free-electron states with energies E < —2c 2 or E > 0, the explicit expressions of the 
radial components Pek(t") and Qek{i") are quite elaborate when compared with the bound 
solutions in Eqs. (|13|) and (|14j) and are not displayed here. They are shown and discussed, for 
example, in the text by Eichler and Meyerhof 0. However, these continuum states are also 
implemented in the Dm AC package and can be accessed as well by means of the procedures 
Dirac_spinor() and Dirac_radial_component(). 



2.1.3 Semirelativistic framework 

The two radial functions (|T3|) and (fTl]) represent, together with Eq. (|1U[). an exact solution 
of the Dirac equation for a pure Coulomb potential. In practice, however, the use of these 
radial components often leads to rather lengthy and complicated computations. Instead of the 
exact form, therefore, the approximate Darwin wavefunctions |19j have been widely used in 
applications such as the theoretical description of the Coulomb ionization 20 or the electron- 
positron pair creation |2j. For a Coulomb potential, —Z/r, the Darwin wavefunctions 

*%L(r) = (l-ia.v)««Wr) (19) 
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are expressed in terms of the nonrelativistic solutions (j2HSJ of the Schrodinger equation 
from subsection 2.1.1, where = (1,0, 0,0) + and u(> = (0, 1,0, 0) + are the two four- 
component spinors for a particle at rest with either spin-up or spin-down, respectively. Ob- 
viously, the Darwin functions Q19JI have a much simpler structure when compared with the 
exact bound-state solutions Q13j) and (|14[) from above; they are accurate to first-order in aZ 
in the expansion of the Dirac equation and are normalized to the same order. 

Apart from the bound states (|19|). approximate continuum wavefunctions can be derived 
which are accurate to aZ in the relativistic motion of the electron and which are known 
as Sommerfeld-Maue wavefunctions [22] from the literature. Again, these functions can be 
expressed in terms of the nonrelativistic functions and have been widely used in applications 
|211 \2'6\ . Although the explicit form of these functions is slightly more complicated than in the 
bound case, their derivation follows very similar lines. Here, we will not display these functions 
in detail but simply refer to the procedures Dirac_Darwin() and Dirac_Sommerfeld_Maue() of 
the Dirac code which facilitate the access to both, the Darwin and the Sommerfeld-Maue 
semirelativistic approximations. 



2.2 Green's function 

Apart from the wavefunctions, which describe the electron in some particular quantum state, 
one often needs the summation over all the (unoccupied) states, especially if parts of the atomic 
interaction are treated as a perturbation. Such a summation over the complete spectrum is 
required, for example, for describing the two- (and multi-) photon processes in hydrogen-like 
ions 24 or for the computation of hyperpolarizabilities. Although, in principle, it appears 
straightforward to carry out such a summation explicitly by using the function from above, 
the large number of terms and the need for integrating over the electron continuum may 
hamper such an approach. Instead, the use of the Green's function 

G{T y ;E) _ ^ H« (20) 

is usually much more favorable since it provides a simple access to the perturbative treatment 
of various types of interactions in the computation of atomic properties. 

As easily shown, the Green's function (|2Uj) can be obtained generally as solution of the inho- 
mogeneous differential equation 



H-E 



G(r,r';E) = 6(r - r') , (21) 



where H denotes the Hamiltonian operator of the system. Similar to the case of the wave- 
functions, the form of the Hamiltonian H and, hence, the explicit expressions for the Green's 
function depends of course on the framework in which the system is described. For the 
Schrodinger operator H = — V 2 /2 — V(r) (the nonrelativistic case), for example, the solu- 
tions of Eq. (|2l"1) can be separated - again - into radial and angular functions by means of 
the ansatz 

1 oo 

G(v,r';E) = —Y J Y.9l(ry-,E)Y^'^')Y lm {'d^). (22) 

1=0 m 
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For a Coulomb-potential V(r) = —Z/r in the Hamiltonian above, the radial functions are 
known analytically in terms of (special) Whittaker functions of the first and second kind 



9i(r, r'; 



E) 



v r(i + 1 - v) 

~Z T(2Z + 2) 




(1r < Z\ f2r > Z\ , , 



where v 



z 



and r> = max(r, r' 



■') and where r. 



min(r, r') refer to the larger and 



smaller values of the two radial coordinates, respectively. 

A similar radial-angular representation can be found for the (relativistic) Dirac Hamiltonian 
in Eq. (|21j) but leads to a quite complicated form for the radial part which is given then by 
a 2 x 2 matrix. In the Dirac program, we provide the procedure Dirac_Greens_radial() which 
supports the computation of the radial functions (|23|) or of the corresponding relativistic 
Green's functions. 

2.3 Integrals with hydrogenic wavefunctions 

The Coulomb wave and Green's functions from the previous subsections can be applied for 
studying the properties and dynamics of hydrogen-like ions. Usually, they are the basic 
ingredients in defining the transition amplitudes and cross sections and, by using their radial- 
integral decompositions in Eqs. @, (|1U|). and 1|22|) immediately lead to a number of radial 
integrals including various types of (transition) operators. While, of course, these integrals 
can always be determined numerically, an analytic (pre-) evaluation of these integrals is in 
many cases possible and results in a considerable gain in the computation of the transition 
amplitudes. Often, moreover, the analytic evaluation of the radial integrals is faster and more 
accurate when compared with a (straightforward) numerical integration. For several integrals, 
which frequently arise in applications, therefore, the analytic formulas have been implemented 
explicitly into the Dirac code. Below, we first discuss the matrix elements containing the 
operators r k , e~ UT , d m /dr m , the integrals for studying the light -matter interaction or for 
collisions with external particles, but turn later also to the (well-known) Slater integrals for 
treating the interaction among the electrons. 

2.3.1 Radial integrals of r— dependent operators 

Nonrelativistic framework: For calculating the structure of the hydrogen-like ions and their 
cross sections in collisions with atoms and electrons, one often needs the matrix elements 
for the operators r k , d m /dr m , e~ UT , or various combinations of these operators. Within 
the nonrelativistic Schrodinger theory, the analytic calculation of such r-dependent matrix 
elements have been analyzed and discussed especially by Sanchez and Lopez Pineiro [25] who 
applied the hypervirial theorem along with a second quantization formalism. They showed 
that these matrix elements can be typically reduced to radial integrals of the form 
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(nl\r k e~ ar \n> ' V) = / P n i{r)r k e~ ar P n >v(r) dr 



o 

^ — ' * — ' to! to'! \n — I — 1 — m / \n — t — 1 — m 

m=0 m'=0 x 7 x 

2Z\ 1 ' / 2Z\ 1 ^" 1 { Z Z ^ — (fe+i+i'+m+m'+3) 



x U + - + - 



x r(fc + l + /' + to + to' + 3) , (24) 

which can then be utilized also (if necessary) in order to compute the matrix elements including 
some derivatives d m /dr m P n i(r) of the wave functions. For the sake of brevity, however, 
we will not display these recursion formulas but only mention that, for all the operators 
from above, the corresponding radial integrals can be evaluated by means of the procedure 
Dirac_r_matrix_element() of the Dirac program. In this procedure, an analytic evaluation of 
the integrals is carried out whenever possible. 

Relativistic framework: Of course, an analytical expression of the matrix elements with oper- 
ators r k e~ ar can be found also for the relativistic Dirac wavefunctions. However, in order to 
derive such expressions it is more convenient not to start from the standard representation of 
the hydrogenic functions ()13H14|) but first to re-write them in terms of a series expansion in r 
as it was suggested originally by Rose PS] . By using such r-series representation, the r k e~ ar 
matrix elements can be written, similar to (|24j). in term of Y- functions as 



(nK|r fc e-° r |n'«') = / (P nl {r) P n , K ,(r) + Q nl (r) Q n , K ,(r)) r k e~° r dr 

o 

n— \k\ n'— \k'\ 

= N nKi N n i K i ^2 X> { C nK,,v C n l k' ,v' + C nn,v c n ' k' ,v<) 



u=0 u'=0 
s+v-l Js'+v'-l 



(q + q' + -) 1 +»+»'+ s + s '+ k 



Tiv + v' + k + s + s' + l), (25) 



where the parameters q, s and the normalization factor N nK were defined by Rose [26], and 
where the coefficients c£ K v are given by 



1/2 ((_„ + !«!)„ 2 ") 



liV^^J " npTT\ [(y-n+\K\)±{aZ/q- K )] (26) 

/ V\ (2S + \)y 



with [a) v being the Pochhammer symbol. 



In contrast to the radial integral (|25j). no analytic expressions for the matrix elements con- 
taining the derivative d m /dr m within the relativistic framework are implemented explicitly 
in Dirac program. For all such matrix elements, however, the (straightforward) numerical 
integration is applied by making use of the standard Maple command int(). 
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2.3.2 Grant's integrals 

Relativistic framework: As known from relativistic collision theory [3JI27], the interaction of 
(atomic) electrons with the radiation field or with other atoms and ions can usually be traced 
back to the computation of the following three radial integrals 

/•oo 

l£(q;a,b) = / drj L (qr) (P a (r)Q b (r) ±Q a (r)P b (r)) (27) 
J o 

/•oo 

l£{q;a,b) = / drj L {qr)P a (r)Q b {r) (28) 
J o 

/•oo 

J L (q;a,b) = / drj L (qr) (P a (r)P b (r) + Q a (r)Q b (r)) (29) 
J o 

where jiiq^) denotes a spherical Bessel function of order L and, as before, P{r) and Q{r) are 
the radial components of the one-electron wavefunctions (|l(Jj) . Because of the oscillations in 
the Bessel as well as the radial wavefunctions, however, an accurate value of these integrals is 
not always easy to obtain numerically. In the Dm AC program, we therefore implemented the 
analytic formulas for these integrals in terms of a finite series over the confluent hypergeometric 
functions as known for the bound-bound and bound-free transitions [281 I29j . While this 
typically results in a much faster (and more accurate) computation of the integrals, it requires 
however that at least one set of radial components — either P a (r), Q a (r) or P b (r), Q b {r) - 
has to represent some bound-state electron. 

Nonrelativistic framework: Since, in the limit c — > oo, the small component Q(r) becomes neg- 
ligible, only contributions from the last integral (|29|) survives in the nonrelativistic description 
of interaction processes with the radiation field and gives rise to the integral 

POO 

J L (q; a,b)= / drj L (qr) P a (r)P b (r) . (30) 
Jo 

An analytic evaluation of this integral with nonrelativistic radial orbital functions is supported 
by the procedure D i ra c_l J _ra d i a I ( ) of the Dirac program for both, bound-bound and bound-free 
electronic transitions. 

2.3.3 Slater integrals 

Although the main emphasis in developing the Dirac program has been placed on studying 
the properties and behaviors of the hydrogen-like ions, the wave functions from above can be 
used also to estimate the effects of the electron-electron interaction for few-electron atoms 
and ions. For such systems, one often needs to compute the (two-electron) matrix elements 
for the Coulomb repulsion operator I/Y12 which, after some algebra, leads to the so-called 
Slater integrals 

/•oo /-co 

R k (a,b,c,d)= dr x dr 2 P nala {ri)P nb i b {r 2 )^P nclc {n)Pn d l b {r2), (31) 
jo Jo r> 

where r> = max(ri, r 2 ) and r< = min(ri, r 2 ) refer to the larger and smaller values of the two 
radial coordinates, respectively. The analytic evaluation of this integral in terms of confluent 
hypergeometric functions has been implemented in the Dirac package for the nonrelativistic 
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framework. For reasons of efficiency, however, no attempt has been made to incorporate 
also the 'relativistic version' of the Slater integrals which would raise much larger demands 
regarding the accurate computation of these values. 

3 Program organization 

3.1 Overview about the DlRAC package 

The Dirac package has been developed as an interactive tool within the framework of Maple 
(version 8). It is designed in order to provide a simple and consistent access to the wave 
functions, the excitation and decay properties as well as the dynamical behaviour of the 
hydrogen-like ions. To allow for a rather large flexibility in the use of the program, most 
of the functions and quantities from the hydrogen atom model are supported both, within 
the nonrelativistic as well as the relativistic theory. With the present version of Dirac, we 
start with compiling different (analytical) representations of the Coulomb wave and Green's 
functions using polar coordinates. In addition, both the analytical and numerical evaluation 
of the - most frequently - required integrals are supported and now facilitate their accurate 
computation for a large number of properties. 

Following Maple's philosophy, the Dirac program is organized in a hierarchical order. In 
the current version, there are about 50 procedures at quite different level of complexity from 
which most however are kept hidden to the user. As in Maple's own design, these procedures 
mainly serve as language elements in order to build up commands at some higher level of the 
hierarchy. Therefore, only about 15 procedures at top level need to be known by the user 
and are briefly summarized in Table 1. More detailed information about the use of these 
procedures, their optional arguments, or the output can be obtained from the user manual 
Dirac-commands.pdf which is distributed together with code. In this manual, we mainly 
follow the style of the help pages of Maple and The Maple Handbook by Redfern [30 from 
earlier years. To illustrate some of the basic features of the program, moreover, a few selected 
examples are displayed below. 

As known from Maple's recent up-grades, most of its internal commands make use of rather 
short names and often of some abbreviations for the various quantities and methods. Although 
such a name convention might be favorable for the frequent use of the commands, these names 
are usually not easy to remember and to recognize in their output if implemented in a code. In 
the Dirac program, therefore, we follow a slightly different concept by making use of rather 
long names to better explain the purpose of the procedures. In the long run, hopefully, this 
will simplify the readability and the maintenance of the program. Moreover, all the commands 
of the Dirac package begin with the additional prefix Dirac_ in order to distinguish them 
from Maple's internal functionality. 

3.2 Getting started with DlRAC 

The Dirac package is distributed as a compressed tar-file dirac. tar from which the dirac 
root directory is obtained by the command tar -xvf dirac. tar. This root contains the 
source code library (for Maple 8), a Read. me for the installation of the package as well as 
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Table 1: Main commands of the Dirac package. A more detailed description of these and a few less 
important procedures is given in the user manual Dirac-commands.pdf and is distributed together 
with the code. 



Procedure 



Description 



Dirac_convert() 
Dirac_energy() 

Dirac_evalf() 

Dirac_Greens_radial() 
Dirac_IJ_radial() 
Dirac_orbital() 
Dirac_r_matrix_element() 

Dirac_radiaLcomponent() 

Dirac_radiaLorbital() 

Dirac_settings() 

Dirac_Slater_radial() 

Dirac_spin_orbital() 

Dirac_spinor() 

Dirac_Ylm() 



Converts the physical results between different unit systems. 

Returns the nonrelativistic energy E n from Eq. (@J or the relativistic 
energy E nK from Eq. (|12|l for bound-electron states. 

Attempts to simplify an algebraic expression for either the energies, 
wavefunctions, or matrix elements by substituting the numerical values 
for some fundamental constants. 

Calculates the radial part of either the nonrelativistic Coulomb Green's 
functions (|23|l or of the corresponding relativistic Green's functions. 

Evaluates the radial integrals (|27Hr?()|) for both, the bound-bound and 
bound-free electron transitions. 

Returns the nonrelativistic hydrogenic orbitals tp n (E) im( r , ^> if) from Eq. 
P|l as function of the given arguments. 

Evaluates the radial matrix elements for various r-dependent operators 
between hydrogenic wavefunctions, both within the nonrelativistic and 
the relativistic framework. 

Returns the large and small radial components, P n (E)n( r ) & n d 
Qn(E) n( r ), of the relativistic orbital functions I(13H14(1 . These compo- 
nents are supported for both, the bound and free-electron states. 



Returns a nonrelativistic radial orbital P r 



(E)l; 



i(r) for both, the bound 



and free-electron states. 

Specifies the framework, nuclear charge and the units which are to be 
used by the DiRAC program. 

Evaluates the Slater integral 1)31(1 with nonrelativistic orbital functions. 
Returns the Dirac spherical spinor f2 Km from Eq. 

Returns the relativistic hydrogenic spinor ip n (E) smf^iJ, ^) from Eq. 

JEJ. 

Return the spherical harmonic Yi m (i), ip). 
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the program manual dirac-commands.pdf. This latter document explains all the user rele- 
vant commands along with the output format, their optional arguments and some additional 
information which might be of interest for the application of the procedures. The dirac root 
directory also contains an example of a .mapleinit file which can be easily modified and 
incorporated into the user's home. Making use of such a .mapleinit file, the module dirac 
should then be available like any other module of Maple simply by typing 

> with(dirac) ; 

Welcome to Dirac 

Dirac_save_nuclear_charge := Z 
Dirac_save_f ramework := nonrelativistic 
Dirac_save_units := atomic 

As seen from the printout above, the settings of the Dirac program are specified by the 
three global variables which must not be overwritten interactively. The (global) variable 
Dirac_save_f ramework, for instance, specifies the theoretical framework in which the var- 
ious commands will respond to the user; it may take the two values relativistic and 
nonrelativistic and determines of how the input and output of the procedures is to be 
treated internally. In practice, of course, the input/output data of the procedures depend 
not only on the framework but often also on the present choice of the units as specified by 
the variable Dirac_save_units. In the DlRAC program, we presently support atomic Hartree 
(h=m e =e 2 / Aireo = 1), natural (m e =c=e=l), SI, and Gaussian units which are recognized by 
setting Dirac_save_units to one of its allowed values: atomic, natural, SI, or cgs, respec- 
tively. Moreover, most of the solutions and output of interest will depend on the nuclear charge 
Z as internally specified by the global variable Dirac_save_nuclear .charge. To change these 
particular specifications of the Dirac program, the command Dirac_settings() can be used; if 
instead of 'atomic' units, the user wishes to employ 'natural' units for the input and output 
of the procedures, the initial specification can be modified by entering 

> Dirac_settings (units , natural); 

Units are changed to natural 

Therefore, by specifying the global variables properly due to the demands of the user, the 
procedures of the Dirac program can be adapted to the particular task. A proper choice of 
these global definitions makes the program appropriate for a wide range of applications. 
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4 Interactive work with the Dirac program 

To illustrate the use of the Dirac program and to provide some test cases, below we will 
display and explain a few examples concerning the fine-structure of the hydrogen-like ions 
as well as the (algebraic) access to their wave functions. A more advanced example later 
shows the computation of matrix elements within the relativistic framework as they arise 
frequently in the study of light-matter interactions. For all these examples, we assume the 
Dirac program to be loaded into the current Maple session as explained in the previous 
section 13.21 

4.1 Fine structure of hydrogen— like ions 

In the nonrelativistic Schrodinger theory, the level energies (j3J) only depend on the principal 
quantum number n and, hence, are degenerate with respect to the orbital quantum number I 
and the magnetic quantum number m. This 2 n 2 degeneracy (the factor 2 arises from the spin 
of the electrons), however, is partially removed if Dirac's theory is applied. In the relativistic 
case (|12j) . namely, the nonrelativistic levels with energy E n split into n different sublevels 
according to the total angular momentum j = 1/2, 3/2, ... n — 1/2 of the state. This splitting, 
which is known as the fine-structure of the hydrogen-like ions, can be easily illustrated below 
by means of the Dirac program. 

To obtain insight into this splitting as a function of the nuclear charge Z, let us first compare 
the relativistic and nonrelativistic level energies for n = 2. Since the nonrelativistic framework 
is predefined whenever the Dirac program is (re-) loaded, we start with the Schrodinger 
energy (J3J) which is returned for an arbitrary nuclear charge Z = Zq by entering 

> Dirac_settings (charge , ZO) ; 

Nuclear charge is changed to ZO 

> E_2(Z0) := Dirac_energy(2) ; 

2 

ZO 

E_2(Z0) := 

8 

Of course, the well-known nonrelativistic energy expression E n= 2 = —Zq/8 is returned, 
given in atomic units as the default setting. Analogue, we may obtain the n = 2 energy 
levels as predicted by the Dirac theory. To this end, we first need to change to the relativistic 
framework for all subsequent calls of the Dirac procedures 

> Dirac_settings (framework, relativistic); 

Framework is changed to relativistic 

and, again, make use of the same command Dirac_energy () . In the relativistic framework, 
this procedure now expects the two arguments n and k to return the level energy E nK from 
Eq. (J12j) . To get the energy of a 2s 1/ / 2 electron with quantum numbers n=2 and k = -1, we 
type 
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> E_2s(Z0) := Dirac_energy(2, -1) ; 



1 

1 

I 2 2 \l/2 

I ZO Dirac_Constant_alpha I 
U + 1 

I 2 2 1/2 2| 

\ (1 + (1 - Dirac_Constant_alpha ZO ) ) / 

E_2s(Z0) := 

2 

Dirac_Constant_alpha 

which returns this energy in the most general form, including the nuclear charge ZO and the 
(unevaluated) fine structure constant a as yet given by the global variable 
Dirac_Constant_alpha. At the first glance, such a general form may look quite complicated 
but it helps of course in treating the expressions algebraically. For our discussion of the 
fine-structure behaviour as a function of Z, we may first assume the aZ parameter small 
when compared with unity and apply the standard Maple command series () to expand 
this expression for the 2s]/ 2 energy in powers of Z 

> series (E_2s (ZO) , ZO) : 

> E_2s_exp(Z0) := convert (%, polynom) ; 

2 2 4 

E_2s_exp(Z0) = -1/8 ZO - 5/128 Dirac_Constant_alpha ZO 

Subtracting from this expression the Schrodinger energy E n= 2 from above 

> E_2s_exp(Z0) - E_2(Z0); 

2 4 

5 Dirac_Constant_alpha ZO 



128 

we find the so-called relativistic Pauli correction AE nK to the (nonrelativistic) level energy 
^ |2] which, as seen from the printout, scales like a 2 Zq and, thus, will strongly increase for 
high-Z ions. 

The result above for the Pauli correction AE nK displays the analytic behaviour of the fine- 
structure splitting; to get a better impression on the size of this splitting, we may evaluate 
this expression for n = 2 also numerically. To substitute the value a = 1/137.036 of the fine- 
structure constant for Dirac_Constant_alpha, we may use the command Dirac_evalf () in 
the Dirac program and obtain 

> E_2s(Z0) := Dirac_evalf (E_2s(Z0)) ; 
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18778 . 865231694104538 

E_2s(Z0) := 

/ 2 \l/2 

I 0.000053251353990881516089 Z0 | 
U + 1 

I 2 1/2 21 

\ (1. + (1. - 0.000053251353990881516089 ZO ) ) / 



- 18778.865231694104538 

By evaluating in a similar way also the level energy for the 2p 3 / 2 shell (n = 2, k = -2) 

> E_2p(Z0) := Dirac_energy(2, -2); 

> E_2p(Z0) := Dirac_evalf (E_2p(Z0)) : 

we can finally calculate the level shifts of the j = 1/2 and j = 3/2 fine-structure levels with 
respect to the value from the nonrelativistic theory. For atomic hydrogen (Zq = 1), this gives 
rise to 

> subs(Z0 = 1, E_2s(Z0) - E_2(Z0)); 

> subs(Z0 = 1, E_2p(Z0) - E_2(Z0)); 

-5 

-0.208018917100000 10 

-5 

-0.41602897100000 10 

in atomic units and can be compared immediately to the values of 0.46 cm -1 (0.208 x 10~ 5 
au) and 0.09 cm -1 (0.416 x 10 -6 au), respectively, as known from the literature 012]. 

4.2 Relativistic contraction of the s— orbitals 

In section 2, we have displayed the radial components for a bound-state electron for both, 
the nonrelativistic © and relativistic theory (|13H14|) , The spatial behaviour of these one- 
electron solutions for different nuclear charges Z has been discussed in detail at a number of 
places to illustrate the 'relativistic effects' for high-Z ions. While, for small charges Z, 
the radial Schrodinger function (J5J) and the large component ()13j) from the Dirac theory are 
practically the same, a comparison of these functions for high-Z ions shows a clear contraction 
of the large component towards the nucleus. This effects, which is known as the relativistic 
contraction of the wavefunctions, influences not only the structure of heavy (hydrogen-like) 
ions but, for instance, also their ionization behaviour ^ST or the chemical binding in materials 
containing heavy elements. 

To illustrate the contraction of the 3s wavefunction for a hydrogen-like uranium ion U 91+ 
by means of the Dirac program, we shall extract and compare the radial function (JSJ) with 
the large component (|13|) from Dirac's theory. The algebraic expression for the P n= z, i=o (r) 
nonrelativistic function can be obtained by 
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> Dirac_settings (framework, nonrelativistic) ; 

Framework is changed to nonrelativistic 

> Dirac_settings (charge , 92); 

Nuclear charge is changed to 92 

> P_3s(r) := Dirac_radial_orbital(3, 0, r) ; 

368 1/2 1/2 92 r 2 

P_3s(r) := — 23 3 exp( ) (27 - 1656 r + 16928 r ) r 

243 3 

where, in the first two lines of printout, we set the nonrelativistic framework and define the 
nuclear charge to be Z = 92. 

Returning to the relativistic framework as before, the large component of a 3s electron is 
generated in the Dirac program by calling 

> P_3s_r(r) := Dirac_radial_component (large, 3, -1, r) : 

for which, we suppress here the printout to screen (using the colon at the end of the line) 
as this would result in quite a lengthy expression in terms of the hypergeometric function 
and various physical constants [cf. Eq. (|13|)]. For the present purpose, instead, we first better 
'evaluate' all the physical constants to their numerical values by calling 

> P_3s_r(r) := Dirac_evalf (P_rel (r) ) ; 

. 741 13463002424792093 

P_3s_r(r) := 69.235282776899202095 r 

exp(-32. 599243267231110742 r) (3.8221513992160292467 
hypergeom([-2.] , [2.4822692600484958419], 65 . 198486534462221484 r) - 

2. hypergeom([-l.] , [2.4822692600484958419], 65 . 198486534462221484 r) ) 

Having the 3s radial wave functions from the nonrelativistic and relativistic theory, we can 
plot them by using the standard Maple procedures plot () . Physically more important than 
the wave function, however, is the 'contraction' in the charge distribution as given by the 
radial distribution functions \P n =3, z=o( r )| 2 an d \P n =3, k=-i(^)\ 2 , respectively. 

> plot([(P_3s(r))~2, abs(P_3s_r(r))~2] , r=0..0.3, linestyle= [DASH, SOLID]); 
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As expected, a clear contraction of the relativistic function (solid line) is found when compared 
with the nonrelativistic radial wavefunction (dashed line). 

4.3 Evaluation of the relativistic radial integrals 

Until now, we have demonstrated the use of the Dirac program by simple 'textbook examples' 
such as the computation of the level energies or the hydrogenic wave functions. The Dirac 
procedures, however, are able to help the user with much advanced studies concerning the 
excitation and decay properties or the dynamical behaviour of the hydrogen-like ions. For such 
investigations, we need to compute transition matrix elements including the wavefunctions 
from above and the corresponding (transition) operators. The exact relativistic treatment of 
the photoionization of hydrogen-like ions, for instance, can be traced back to the computation 
of the bound-free transition amplitude [3] 



and the final continuum electron, respectively, and where the operator u\e lkr describes the 
electron-photon interaction. 

The simplification of the transition amplitude (|H2|) requires elements from the algebra of 
angular momentum and has been discussed elsewhere jSlE]- I n particular, it was shown that 
this amplitude can be factorized into radial and angular parts where the latter ones evaluate 
to Wigner 3 — j and 6 — j symbols and a constant factor, while the radial part is reduced 
to the integrals (|27H29j) . Below, we therefore discuss how the Dirac package can be utilized 
to compute these radial integrals. We start with the straightforward numerical integration of 
the radial integral Ii = i(q',a,b) , which arises if a lS]/ 2 electron (n^ = 1, « 6 = -1) is ionized 
by a photon with energy = 70 eV into a Sx/ 2 continuum state (partial s— wave; Kf = -1). 
For this particular choice of the parameters and quantum numbers, we may simply enter 

> Dirac_settings (charge , 1): 

> Dirac_settings (framework, relativistic): 




(32) 




(relativistic) wavefunctions for the initially bound 
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> Pa(r) := Dirac_evalf (Diracjradial_component (large, 1, -1, r)): 

> Qb(r) := Dirac_evalf (Dirac_radial_component (small, 2.072432345, -1, r,/ree)): 

> int(Pa(r)*Dirac_spherical_j (1 , . 01877198124*r) *Qb(r) ,r=0. .70. ) ; 

-5 

-0.23633058961153019500 10 

where the keyword free in the procedure Dirac_radial_component() is used to generate the small 
component of a continuum state with (kinetic) energy Et = 2.072432345 (au) and where most 
of the intermediate printout has been suppressed by using a colon to terminate the first 5 lines. 
Both the energy Ef and the momentum transfer q = 0.01877198124 (au) are easily derived 
from the conservation of energy, assuming / = 0.5 (au) as the ionization potential for the Is 
electron of atomic hydrogen. In the last input line, moreover, the procedure Dirac_spherical_j() 
is used to returns the spherical Bessel function jL^or). 

The (direct) numerical computation of the integrals (|27H29j) is typically reliable but also very 
time consuming. For the integral above, for instance, Maple requires about 145 seconds at 
a 1 GHz processor. An alternative and (much) faster computation is achieved by using the 
analytical expressions for the radial integrals (|27H29|) in terms of the special hypergeomet- 
ric functions [2Hj- Within the DlRAC program, such a (analytical) integration algorithm is 
supported by means of the procedure Dirac_IJ_radial() which can be invoked with the same 
parameters as above to determine the /£ =1 (g;a,6) integral 

> Dirac_IJ_radial(IL0, 0.01877198124, 1, 1, -1, 2.072432345, -1, free): 

-5 

-0.23633058961153083452 10 

Although, of course, the integral is the same for using both, the numerical or analytical 
integration scheme (at least up to the 16-th decimal digit), the procedure D i ra c_l J _ra d i a I ( ) is 
faster by almost two orders of magnitude when compared to Maple's numerical integration. 
It is also more accurate as found by increasing the number of digits to Digits := 30. 

5 Summary and outlook 

The Dirac program has been presented to facilitate a fast and consistent access to the basic 
elements and properties of the hydrogen-like ions. In its present form, this program consists 
out of a set of Maple procedures for the numerical and analytical evaluation of the Coulomb 
wave and Green's functions, both within the nonrelativistic and relativistic framework. Apart 
from these functions, moreover, the computation of radial integrals (as they frequently arise 
in the evaluation of transition amplitudes for various types operators) is supported by means 
of the Dirac program by implementing the analytic formulas. 

Of course, the procedures above can be seen also as the technical basis for more advanced 
studies concerning the behaviour of the hydrogen-like ions or, more generally, applications of 
the atomic shell model. In the future, therefore, several lines are possible (and desirable) for 
the further development of the Dirac program. In a first step, we want to extent the present 
set of procedures towards the interaction of hydrogen-like ions with the radiation filed as it 
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occurs in the ionization or capture of electrons by high-Z ions or in laser fields. Examples of 
very desirable procedures refer to a simple access to the decay rates, angular distributions, or 
to the polarization of the emitted radiation. Having collaborations with the atomic physics 
community at GSI and elsewhere for many years, we know that these developments will make 
the Dirac package attractive to both, experimental and theoretical investigations on energetic 
ion-atom and ion-electron collisions. 
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